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Abstract 

We study the stability of the 0(N) fixed point in three dimensions under 
perturbations of the cubic type. We address this problem in the three cases 
N = 2,3,4 by using finite size scaling techniques and high precision Monte 
Carlo simulations. It is well know that there is a critical value 2 < N c < 4 
below which the O(N) fixed point is stable and above which the cubic fixed 
point becomes the stable one. While we cannot exclude that N c < 3, as 
recently claimed by Kleinert and collaborators, our analysis strongly suggests 
that N c coincides with 3. 
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1 Introduction 



Quantum field theories with 4 type interactions are of importance in several phys- 
ical contexts. In particular, they represent one of the most powerful tools in the 
study of critical phenomena U . Due to their simplicity they allow perturbative ex- 
pansions up to rather large orders from which one can extract estimates for various 
critical quantities (critical indices and amplitude ratios) comparable in precision 
with those of the most advanced Monte Carlo simulations. In the simplest case 
the theory contains a single field and describes the Ising universality class (for a 
recent comparison between field theoretic and Monte Carlo predictions in this case 
see for instance 0). 

When the field has more than one component the situation becomes more 
complex and different quartic interaction terms can be defined. The simplest one has 
the form (X^i < P 2 ) 2 - It is 0(N) symmetric and describes the O(N) universality class 
to which belongs, for instance, the isotropic ^-component Heisenberg ferromagnet. 
Besides this term, the most interesting additional contribution is YliLi (fit which 
breaks the O(N) symmetry but preserves the cubic invariance. The cubic subgroup 
of O(N) is composed of permutations and reflections of the N components of the 
field. Note that in the following "cubic" always refers to the symmetry and not to 
a third power. The importance of the cubic term is due to the fact that in a real 
crystal the crystalline structure gives rise to anisotropies which are mainly of the 
cubic type. Thus real crystals are better described by mixed actions in which both 
the O(N) and the cubic term are present. 

Besides this phenomenological reason, this mixed model is also interesting in 
itself as it is a simple non-trivial QFT with different fixed points in competition 
among them. In fact, it is easy to see that in this model there are four possible fixed 
points: the trivial Gaussian one, the Ising one (which corresponds to the situation 
in which the N components 0, decouple), the O(N) symmetric and the cubic one 
(see fig.l). It was shown more than twenty years ago || that while the Gaussian 
and Ising fixed points are always unstable, the O(N) and cubic ones interchange 
their role as N increases. For N < N c the O(N) symmetric point is stable, while 
for N > N c it is destabilized by the cubic interaction and the cubic fixed point 
becomes the stable one (see fig.l). It is possible to see within the framework of 
the e-expansion that in three dimensions N c < 4. The common lore, (supported by 
e-expansion up to the third order) has always been that N c should lie somewhere 
between 3 and 4 in three dimensions, thus implying that the N = 3 case, which is 
the most interesting one for applications to real crystals, should have a stable 0(3) 
symmetric fixed point. 

In these last years, this commonly accepted scenario has been contrasted in a 
series of papers ||, ^J, [| [7J Q which suggested that N c should lie below 3. As a 
consequence the critical behaviour of magnetic transitions in real crystals should be 
described by the cubic symmetric fixed point, a result which, if confirmed, would 
be of relevant interest from a theoretical point of view. 
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The aim of this paper is to test this conjecture with a high precision Monte Carlo 
simulation. By studying finite size corrections of a cubic invariant perturbation term 
exactly at the critical O(N) point we can extract the eigenvalues of the stability 
matrix of the O(N) fixed point. We study the three interesting cases N = 2,3,4. 
For N = 2 and N = 4 the expected results (stability of the isotropic and cubic 
fixed point respectively) are immediately visible from the data. In the N = 3 
case our results imply that N c m 3. Obviously the numerical simulation can not 
decide whether iV c = 3 is an exact result. However we obtain an upper bound for 
the absolute value of the stability index |&2| for the 0(3) fixed point, which turns 
out to be impressively small. In particular we are able to exclude all the existing 
estimates (|, |5j, || [K| [11], [TJJ except that of Kleinert and collaborators || |7|, ||, 
which is still compatible, within one standard deviation, with our result. 



2 The cubic model 

We are interested in the three dimensional quantum field theory defined by the 
Lagrangian 

1 N \ N ,1 N 

z »=1 4 - i=l 4 - i=l 

While the term (5Zi=i 0?) 2 is 0(N) symmetric, the term Y,i=i 4>\ is only invariant 
under the "cubic" subgroup composed by permutations and reflections of the N 
components 4>i- This model is discussed in great detail in several quantum field 
theory text books |]. 

In the following we shall review the first order results in the e-expansion. This 
rather simple approximation already gives all the qualitative features of the renor- 
malization flows of the model. 

The fixed points of the theory are given by the zeros of the /5-functions. The 
stability matrices are given by the derivatives of the /3-functions at the zeros. From 
the eigenvalues of these matrices it is then easy to identify the stable fixed point. 

The two /^-functions are given, at the first order in the e-expansion, by 

2^ + 8 

p u = —eu + u — h uv (2) 



3 

(3 V = -ev + -v 2 + 2uv , (3) 

where u and v are the renormalized couplings related to A and [i respectively. 

By looking at the zeros of the /3-functions one sees that there are four possible 
fixed points: 

1] The Gaussian fixed point u — 0, v — . 

2] The Ising fixed point u = 0, v — |e . 
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3] The Heisenberg (0(N) invariant) fixed point u 




e, v — . 



4] The Cubic fixed point u = j^e, v = 




The stability matrix is defined as 



B = 



( 



d(3 u (u,v) 

du 
d/3 v (u,v) 

du 



d/3 u (u,v) 

dv 
df3 v (u,v) 

dv 



) 



(4) 



At the first order of the e-expansion one obtains 




-e + 



N+8 
3 



U + V 



U 



) 



(5) 



2v 



-e + 2u + 3v 



The corresponding eigenvalues, evaluated at the four fixed points are: 

1] Gaussian : bi = b 2 = — e . 

2] Ising : b x = -|, b 2 = e . 

3] Heisenberg : bi = e, b 2 = • 

4] Cubic fixed point: b\ = Ne, b 2 = • 

It is easy to see that the Gaussian and Ising fixed points are always unstable, 
independently from the value of N. In particular the Ising f.p. has only one 
direction of instability, while the Gaussian one is unstable in both directions. 

The cubic and O(N) fixed points interchange their role as a function of N. For 
iV smaller than a critical value iV c (which at this order in the e-expansion turns out 
to be 4) the Heisenberg fixed point is the stable one and defines the universality 
class toward which the system flows in the infrared limit. 

For iV > N c , b 2 evaluated at the Heisenberg point becomes negative while b 2 
evaluated at the cubic f.p. becomes positive and the cubic fixed point becomes the 
stable one. 

The renormalization flows corresponding to these two situations are reported in 
fig. 1. For N < N c all initial points with u > and v > ^f^u will flow toward the 
0(iV)-invariant, Heisenberg fixed point. For N > N c all initial points with u > 
and v > will flow in the infrared limit toward the cubic fixed point which, for 
N < N c lies in the v < half plane, exactly at iV = N c crosses the v — axis and 
moves for iV > N c in the v > region. 

Initial points outside the above defined regions flow away toward more negative 
values of u and/or v and finally reach the region in which the positivity condition 
for the quartic potential is no longer satisfied. These trajectories are related (from 
the statistical mechanics point of view) to realizations of the cubic model in which 
the phase transition is of the fluctuation-induced first order type. These models 
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have recently attracted much interest as a laboratory to study arbitrary weak first 



order transitions 13 



The last remaining point is now to find the value of N c in three dimensions. 
It is easy to see by looking to higher orders in the e-expansion, or with the help 
of Monte Carlo simulations, that for N = 2 the Heisenberg fixed point is stable 
and that on the contrary for N = 4 the cubic fixed point is the stable one. Thus 
2 < N c < 4. However it turns out to be hard to decide whether N c is greater 
or lower than 3. Equivalently one can look at the sign of the hi eigenvalue at the 
Heisenberg point for iV = 3. If 62 is positive, then N c must be greater than 3. In the 
past twenty years much efforts have been devoted to settle this question. The first 
result was reported in || where the e-expansion for iV c was extended up to the third 
order leading to the estimate N c = 3.128. In agreement with this estimate (but, 
using a completely different approach), Grover, Kadanoff and Wegner [10] obtained 
62 = 0.053 at iV = 3. Few years later, with different approximation techniques, the 



two contrasting results: N c ~ 2.3 |ll[] an d N c ~ 3.4 |12| were obtained. Ten years 
later in ||, |5| a value N c < 3 was suggested. In particular in g, by means of a 
three loop calculation directly in d = 3, the values N c = 2.91 and h 2 = —0.008 were 
proposed. Finally, more recently, Kleinert and collaborators pushed the e-expansion 
up to the fifth order [[/] and obtained a similar answer. First, in they found, (with 
a [2,2] Pade' approximant) iV c = 2.958. Then in ||, by using a careful resummation 
procedure of the fifth order series, they obtained the value 62 = —0.00214 for the 
stability eigenvalue at iV = 3. Due to the nature of these results it is very difficult 
to add sensible error-bars to these estimates. However, it is clear from the above 
discussion that the existing estimates for N c are scattered around N c = 3 and that 
as the various techniques and approximations become more and more refined the 
corresponding estimates for N c get closer and closer to N c = 3. 



3 The simulation 
3.1 The model 

The cubic model discussed in sect. 2 has a simple and straightforward lattice real- 
ization, defined by the action: 

N 

s = -pY, s.«v-a*££(4)* , (6) 

<xy> x i=l 

where s x is a unit vector in R . < x,y > denotes a pair of nearest neighbour 
sites on the lattice. We consider a three-dimensional cubic lattice of size L and 
lattice spacing a = 1. For /i = we have the standard O(N) invariant (Heisenberg) 
model, while for \i 7^ the cubic- invariant perturbation J2 X J2iLi( s x ) 4 breaks the 
O(N) symmetry. In the following we shall study this model in the three cases 
iV = 2, 3, 4. We shall concentrate our main efforts in the N = 3 case. 
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In three dimensions, for \x = 0, the O(N) model undergoes a second order phase 
transition for some value f3 c (which depends on N) of the coupling. In the vicinity of 
such a point the continuum limit can be taken, leading to the O(N) symmetric QFT 
(corresponding to the v — axis in fig.l) discussed in the previous section. The 
presence of such a continuous phase transition is obviously a mandatory condition 
for the whole analysis. 

The simplest way to extract the value of N c from a lattice simulation is to 
determine the stability eigenvalues 62 for N = 3. From the lattice point of view 
the 62 eigenvalue appears as the critical index which controls the behaviour of any 
cubic-invariant (but O(iV)-violating) observable in the vicinity of the Heisenberg 
transition point. 

The most efficient way to evaluate such a critical index is to look at the finite 
size dependence (as a function of the lattice size L) of a suitable observable (to be 
defined below) evaluated exactly at the critical point (3 C . To this end it is necessary 
to have a very good estimate of the critical coupling. Fortunately, (3 C is known with 
very high precision in each of the three cases N = 2, 3, 4 in which we are interested. 
This is one of the reasons for which we have chosen this particular lattice realization 
of the cubic model. 

In tab. 1,2,3 we have collected the most recent results for (3 C both from Monte 
Carlo simulations and from series expansions. HT-6 1 indicates the biased resumma- 
tion of the HT series in which the value of the index 9 is given as input parameter. 
It is interesting to see that all the estimates agree within the errors. 

In tab. 4 we report the values (chosen from tab. 1,2,3) that we used in our simu- 
lations. 

Table 1: Results for f3 c given in the literature for N = 2 



ref. 


method 


Pc 




14} 

14| 
15] 
161 


MCRG 
MC 
MC 
HT 


0.45420(2) 
0.454170(7) 
0.454165(4) 
0.45419(3) 



Table 2: Results for (5 C given in the literature for N = 3 



ref. 


method 


Pc 


I7| 


MC 


0.6930(1) 


Pfl 


MC 


0.6931(1) 


19] 


MC 


0.693035(37) 


15] 


MC 


0.693002(12) 


16| 


HT 


0.69303(3) 


16] 


HT-# 


0.69305(4) 
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Table 3: Results for (5 C given in the literature for N = 4 



ref. 


method 


Pc 




2C| 


MC 


0.9360(1) 




15 


MC 


0.935861(8) 




16 


HT 


0.93589(6) 




m 


HT-6 


0.93593(6) 



N 


Pc 


2 


0.454165(4) 


3 


0.693002(12) 


4 


0.935861(8) 



Table 4: Values of f3 c used in this article. 



3.2 Observables 

There are four natural observables in the model 
in the action: 



E 



<xy> 



and 



The two terms which appear 

(7) 



X i 

The total magnetization, which is the order parameter of the transition: 



(9) 



and the ratio 



R 



E,(M*) 4 
(M 2 ) 2 



(10) 



which quantifies the violation of the O(N) symmetry in the model []. 

In order to study the stability of the fixed point we are actually interested in 
the derivative of < R > with respect to // at (3 = /3 C and \x = 0: 



Di 



d < R > 



1 Other choices are possible for this last observable. For instance the term 

M 2 

V max 



(12) 



(11) 



where M^ ax is the maximal square of a component of the magnetization would work equally well. 
However it turns out that the ratio R defined above is the one which can be measured in the most 
efficient and simple way. 
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In fact this derivative measures the (generalized) susceptibility of the system 
with respect to a cubic perturbation, just like what happens for the ordinary mag- 
netic susceptibility in the case of a magnetic perturbation or for the specific heat 
in the case of a thermal perturbation. At this point, a standard finite size scal- 
ing analysis tells us that the critical index which measures the infrared stability 
of the system with respect to the above perturbations, also controls the finite size 
behaviour (namely the L dependence) of the corresponding susceptibility exactly at 
the critical point. In particular, in the case in which we are interested, Dr should 
behave for (3 = (3 C and /x = 0, as 

£> fl oc IT 6 " . (13) 

where 62 is exactly the stability eigenvalue of the Heisenberg fixed point that we 
are looking for. 

Scaling laws of the type ( |13"D are expected to hold for sufficiently large values 
of L. For small lattices, correction to scaling terms should be expected. It is thus 
important to have results with small statistical errors for large lattice sizes. Few 
preliminary tests have shown us that (at least for the case N = 3) sizes up to 
L = 32 are needed to extract reliable estimates of 62- This requirement represents 
the major technical problem of this work. 

There are two possible choices to compute Dr. 

• It can be computed directly in the simulation at /x = as 

D R =(PR)-(P)(R) . (14) 



• It can be computed by using the finite difference method, i.e. by simulating 
the model at small non-zero values of /x: 

s AO,) - bho (15) 

2/x 

It is easy to recover the relationship between Dr and Z)r(/x). First, let us notice 
that on a finite lattice R(fi, (3) must be an analytic function of its parameters. This 
holds also at (3 = (3 C . Therefore we can Taylor-expand R(fM, (3) in powers of /x for 
fixed (3 = (3 C . Taking the symmetric difference we obtain 

D r = D r M + ~±ii* + 0(ji a ) . (16) 

Both these definitions have their drawbacks. The /x = simulations are af- 
fected by a strong enhancement of the variance (hence of the statistical errors) as 
L increases. It turns out that the statistical error of D R at a fixed number of mea- 
surements increases roughly as L 3//2 . As a consequence too large samples are needed 
to keep the error sufficiently small for L > 16. 
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On the contrary for Djt(fi) at fixed /i the statistical error does not increase 
with L. However the 0(/x 2 ) corrections do increase with L. Reducing these 0(fi 2 ) 
corrections requires to reduce the value of \i. This in turn requires to increase the 
number of measurements to keep the statistical error fixed. 

In the following section we shall discuss a way out of these problems. By using 
the global O(N) symmetry of the model at fi = an improved version of Dr can be 
constructed. This improvement does not change the L dependence of the variance, 
but gives a significant reduction of its magnitude, thus allowing to reach, with a 
reasonable CPU time, lattices sizes as large as L = 32 which are large enough 
to extract the finite size behaviour with the required precision. Most of the data 
that we shall discuss in the last section have been obtained by using this improved 
observable. We also performed, as a cross check, some simulations at finite //. The 
agreement that we find between the values of Dr obtained in these two ways is a 
non trivial check of the reliability of our results. 

3.3 Variance reduced estimator for Dr. 

Variance reduced estimators have the same expectation value as the corresponding 
standard estimators. However their variance is reduced, which allows for more 
accurate results in Monte Carlo simulations than the standard estimator. A general 
principle to construct variance reduced estimators is to look for degrees of freedom 
which can be integrated out analytically. 

In order to obtain a variance reduced estimator of a< ^j > we integrate P,R and 
PR over the global O(N) rotations. 

This is trivial for P and R, but it is less simple in the case of the PR component. 
Since P is a sum over all lattice sites we can commute the integration over the global 
rotations and the summation over the lattice sites. The integral to be solved is hence 
given by 



where T is an element of O(N), DT the Haar- measure, s x the spin at the site x, 
and m a unit vector in the direction of the global magnetization. For symmetry 
reasons the integral only depends on the angle between s x and m, which we define 
as follows: 



The integral (|i~7l) can be evaluated explicitly for any value of N. Details of the 
calculation are reported in the appendix. Here we only list the results in the three 
cases in which we are interested: 




(17) 



ms x = cos(a) 




N 



2 



I 



9 



— cos (4a) 
32 v ; 



(19) 



16 
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N = 3 



/=^co S (4 tt ) + i ^cos(2 tt ) + lg (20 ) 

N = 4 

/ = — CO s\a)- — cos 2 (a) + — (21) 

25 v ; 50 v ; 200 v ; 

3.4 The Monte Carlo Algorithm 

Due to the different symmetries in the models, we had to use different algorithms 
in the two cases /i = and // ^ 0. 

3.4.1 /i = 



In the case /j = we used the single cluster algorithm of U.Wolff [£l|] and the 
microcanonical overrelaxation algorithm. The basic idea of the cluster-algorithm 
is to construct conditional Ising models. This is achieved by allowing only the 
sign-change of the spin component parallel to an unit vector r in R N . The delete 
probability depends on the pair of lattice sites and is given by 

p d {x,y) = min[l,exp(-2(rs x )(rs v ))] , (22) 

where x and y are nearest neighbour sites on the lattice. The vector r is chosen 
with a probability density uniform on S*^ -1 . For each update a new r is chosen. 

The variance of the improved estimator of Dr is mainly caused by local fluctu- 
ations of the spins. Hence it is useful to supplement the cluster-algorithm with a 
fast local algorithm to produce local changes of the configuration. For that purpose 
we used the microcanonical overrelaxation algorithm. 

The elementary update of the algorithm is given by 



= 2 (n x s) 
nz 



where n x is the sum of the nearest neighbour spins of s x . Since neither a random 
number nor the evaluation of the exponential function is needed for this update 
the CPU-time required is rather small compared with a Metropolis or a heat-bath 
update. 

The whole update cycle that we used in our simulations consists of a mixture of 
9 overrelaxation sweeps and K cluster-updates . K is chosen such that the number 
of spins updated in K cluster-updates is of the order of the number of lattice-sites. 
A measurement is performed after each overrelaxation sweeps and after the cluster- 
updates. There are hence 10 measurements in the 9 overrelaxation sweeps and K 
cluster- up dates cycle. 
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3.4.2 /i^O 



For pi 7^ some modifications are needed. We have restricted the vector r such that 
X^(s 1 ) 4 is not changed by the update. This is guaranteed if the sign of a component 
is changed or two components are exchanged or a combination of both. This means 
that r is either parallel to an axis or is diagonal in a plane, 
r is in 

(1,0,...,0),(0,1,...,0)...(0,0,..,1) (24) 

or 

-^=(1, 1, 0), -J=(l, 0, 1), ...-L(0, ...1, 1) (25) 

or 

-L(l,-1,...,0),-^(1,0,...,-1),...-^(0,...1,-1) . (26) 

While this restriction of r to a discrete subset of S^ -1 does not violate detailed 
balance it means that the cluster-update by itself is not ergodic. 

In order to restore ergodicity we supplement the cluster-update with an (ergodic) 
Metropolis update. For performance reasons we also added a local reflection update 
that is microcanonical for /i — 0. A spin is reflected at the sum of its neighbours. 

4 = - 2s x . (27) 

where S x is the sum of the spins on nearest neighbour sites of x. For /i ^ the 
proposal s' x is accepted with a probability 

P acc = min[l,exp(/i;£[(^) 4 - (4)1)] • (28) 

i 

A whole update cycle consists of one Metropolis update, one local reflection 
update plus K cluster-updates. K is chosen, as in the /x = 0, case such that the 
number of spins updated in K cluster-updates is of the order of the number of 
lattice-sites. A measure is performed in each update cycle. 



3.5 Statistical and systematic errors 

We evaluated statistical errors with the standard binning method. Both in the 
/i = and /i cases bins of 1.000 update cycles are chosen (corresponding to 
10.000 and 1.000 measurements respectively). This binning was already performed 
during the simulation since not all individual measurements could be stored on disc. 
Besides the statistical uncertainty we have to face also the systematic error due to 
the uncertainty A[3 C in the estimate of f3 c . To evaluate this error we also measured 
in the simulation the expectation value < ER >. The difference 

< ER> - < E >< R> (29) 
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gives an estimate of the derivative of < R > with respect to (3, from which we can 
obtain the systematic error induced on R by the uncertainty in f3 c : 

(< ER > - < E >< R >) * Af3 c . (30) 

With a similar construction we obtain the error induced on Dr. We can consider 
this as a lower bound on the accuracy that we can reach for the derivatives. It makes 
no sense to reduce the statistical error of Dr below this bound. This observation 
fixes the typical sample size for the simulations, which turned out to be of the order 
of 20.000 bins. 



4 Results and discussion 

4.1 Results at fi = 

We simulated the models with N = 2, 3, 4 at fi = 0, (3 = /3 C , in the ranges 
L 6 [4 - 16] for iV = 2; L G [4 - 32] for iV = 3 and L <E [4 - 20] for iV = 4. The 
results are reported in tables 5, 6 and 7. In the first column we report the lattice 
size and in the second the values of the derivative Dr. the first error in parenthesis 
denotes the statistical uncertainty, while in the second parenthesis the error induced 
by the uncertainty in f3 c is reported. In the last column we report the sample size 
(number of bins times number of measurements in each bin). 



Table 5: Results for Dr in the N = 2 model 



L 


Dr 


statistics 


4 


.011046(20)(1) 


10000 * 10 4 


6 


.010506(34)(2) 


10000 * 10 4 


8 


.010121(51)(2) 


10000 * 10 4 


10 


.009728(58) (4) 


15000 * 10 4 


12 


.009523(76)(5) 


15000 * 10 4 


16 


.008904(117) (4) 


15000 * 10 4 



4.2 Results at \i ^ 

As a test of the above results we also performed some simulations at /i ^ 0, both 
for N = 2 and N = 3. We evaluated the finite /i estimators Dr(^) by using 
eq.([15D. The results are reported in tab. 8 where, in the last line, we also reported 
for comparison the corresponding fi = estimates. 

The agreement between the results obtained with the two approaches is very 
good and makes us confident on the reliability of the /i = set of data. 
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Table 6: Results for Dr in the N = 3 model 



L 


Dr 


statistics 


4 


.019672(19)(3) 


10000 * 10 4 


6 


.020118(19)(5) 


15000* 10 4 


8 


.020187(20)(7) 


20004 * 10 4 


10 


.020196(25)(10) 


20870 * 10 4 


12 


.020152(32)(13) 


21500* 10 4 


14 


.020178(40)(15) 


20770 * 10 4 


16 


.020233(49)(16) 


20750 * 10 4 


20 


.020094(68) (22) 


20150* 10 4 


24 


.020178(84)(28) 


23145 * 10 4 


32 


.020265(140)(39) 


19560 * 10 4 



Table 7: Results for Dr in the N = 4 model 



L 


Dr 


statistics 


4 


.023474(15)(1) 


10000 * 10 4 


6 


.025287(16)(3) 


10000 * 10 4 


8 


.026294(19)(4) 


10000 * 10 4 


10 


.027097(23) (6) 


10000 * 10 4 


12 


.027742(27) (8) 


10600 * 10 4 


16 


.028774(39) (13) 


10000 * 10 4 


20 


.029605(52)(38) 


10050 * 10 4 



4.3 The b 2 index 

We fitted the data obtained at /i — with the scaling law 

D R = CL- b2 . (31) 

The fit results are collected in tab. 9 where in the second column we give the mini- 
mum value L min of L taken into account in the fit. In the third and fourth column 
we report the reduced x 2 an d the confidence level respectively. Finally, in the last 
two column the best fit values of C and 62 are reported. As usual we give in the 
first parenthesis the statistical error and in the second the error induced by f3 c . The 
various fits are plotted and compared in fig. 2 and 3. 

The large value of x 2 clearly indicates that for any value of iV the sample at 
L = 4 is strongly affected by correction to scaling terms and must be discarded. Fits 
without L = 4 have an acceptable x 2 ■ However this fact does not necessarily imply 
that it is justified to ignore corrections to scaling. Hence we regard the fits with 
L m in = 8 as our final result. Still it remains difficult to quantify the systematic 
error due to corrections to scaling. Based on the experience with the finite size 
scaling analysis of other exponents of the Heisenberg model we expect them to be 
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Table 8: Results for Dr for \x ^ 





N = 2, L = 8 


AT = 2, L = 12 


N = 2, L = 16 


AT = 3, L = 12 


M = 4 
H = 2 
H=l 
[i = 0.5 
H = 0.25 


.010117(28)(1) 
.010112(57) (1) 
.010054(74) (1) 


.010503(140)(3) 
.010008(122)(3) 
.009495(118) (3) 
.009505(103)(3) 


.009139(103)(4) 


.020977(64) (9) 
.020306(130) (9) 
.020076(257)(9) 


n = 


.010121(51)(2) 


.009523(76)(5) 


.008904(117)(4) 


.020152(32)(13) 



of the same order of magnitude as the statistical error of the L min = 8 fits. 



Table 9: Results for C and b 2 



N 


Lmin 




C.L. 


C 


b 2 


2 


4 


2.01 


9% 


0.01335(9)(1) 


0.1362(40)(3) 


2 


6 


1.17 


32% 


0.01381(24)(1) 


0.1519(84)(4) 


2 


8 


0.85 


43% 


0.01445(54)(2) 


0.1711(166)(6) 


3 


4 


26.5 


0% 


0.01936(4)(2) 


-0.0174(10)(6) 


3 


6 


1.32 


23% 


0.02005(6)(2) 


-0.0026(14) (7) 


3 


8 


0.71 


64% 


0.02022(10)(3) 


0.0007(20)(9) 


4 


4 


77.5 


0% 


0.01928(3)(1) 


-0.1473(7)(4) 


4 


6 


1.22 


30% 


0.01997(5)(2) 


-0.1321(10)(5) 


4 


8 


0.51 


67% 


0.02008(8)(3) 


-0.1299(16)(8) 



4.4 Discussion and comparison with other estimates 

As it can be seen from tab. 9, our results for = 3 are certainly incompatible 
with all the existing estimates |Q, [5|, || [10], [TT], |TJ|, except that of Kleinert and 
collaborators |], [7], As a matter of fact, if we keep in the fit for A^ = 3 also the 
L = 6 sample we find an impressive agreement with the result 62 = —0.00214 of 
ref. 0. However, as mentioned above, we strongly suspect that the L = 6 sample is 
still affected by correction to scaling terms and prefer to quote as our best estimate 
the L = [8 - 32] result b 2 = 0.0007(20) (9), which is still compatible with the result 
of ref. , but suggests that N c could indeed exactly coincide with 3. In this respect 
it must also be noticed that the trend of the perturbative estimates of b 2 quoted 
in as a function of the order in the perturbative expansion also suggests that A^ c 
converges to 3 in agreement with our result. 

In any case, let us stress again that it is obviously impossible to decide by 
means of a numerical simulation if A^ c = 3 is an exact result and that the fact 
that the difference \N C — 3| is so small, and compatible with zero, might well be a 
coincidence. However we think that it would be worthwhile to look for an argument 
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which explains why the cubic and Heisenberg fixed point in three dimensions should 
coincide exactly for N=3. 
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Appendix 

In this appendix we evaluate the integral 

i = j o{N) (Et(^)i 4 ) (EK T ™)i 4 ) > (a- 1 ) 

where T is an element of 0(N), DT the Haar- measure, s x the spin at the site x and 
m a unit vector in the direction of the global magnetization. For symmetry reasons 
the integral only depends on the angle between s x and m which is defined as 

ms x = cos(a;) . (A. 2) 

We write the integral as 

I{a) = j 3N _ i 6s{^st)j 3N _^{^) , (A.3) 

where J s n-i denotes the integral over the iV dimensional sphere. J s n-2 jS denotes the 
integral over the N — 1 dimensional subspace defined as the set of all the vectors t 
that for any fixed s satisfy the equation st = cos(o;). We choose the normalizations 
so that J s n-i ds = 1 and J s n-2 s dt — 1 

Because of symmetry we can restrict the calculation to the first component of s 

'<«>-* / s „i»!/ s „/'E';i ■ (a.4) 

Now we decompose the integral J s n-i into the integral over the si component 
and for fixed s\ over the remaining S N ~ 2 . 
We get 

1(a) = N const f ^ d Sl (1 - s?)(^3)/2 g 4 f dg , f df /v-^ ^ 

J Sl =0 Js N - 2 Js N -2,s i 

where 

const, = \ [ Sl=1 dsi (1 - slY N ^ 2 } (A.6) 

Jsi=0 

and s' is s without the 1-component. 
Let us now study 



ds' / dt . (A.7) 

This measure for t is invariant under rotations around the 1-axis. The non-trivial 
question is the measure for the 1-component of t. The range of t\ is given by 



tmax = cos(a)si + sin(aWl - sf (A. 8) 
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and 



tmin = cos(a)si - sin(a)yi - sf . 



(A.9) 

The measure between this extreme values is given by the fact that for any s, t is 
distributed on a S N ~ 2 sphere. Hence the measure is (for N > 2) 



const 2 



t\ — c 
2s 



(JV-4)/2 



(A.10) 



with c = cos(a;)s 1 and s = sin(a)y / l — s\ . The normalization const 2 is given by 

(A.ll) 



const-. 



r-c+s 




(JV-4J/2 ^ 






J t\=c—s 







For fixed ii the integration of the remaining components gives us 
(1 - t\) 2 < R > N _ U with < R > N = 3/(N + 2). 

Now we are in the position to write down the full integral: 



Ka) 



si=l 



'1; 



N consti I d«i(l - -s 

si=0 

(JV-4)/2 



2n(7V-3)/2 4 



'1 



const 2 



c+s 

tl=C — S 



1 - 



2s 



iV+ 1 



(A.12) 
(A.13) 



This integral can be solved with standard techniques and yields in the three 
cases N = 2, 3, 4 the results listed in sect. 3 
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Figure Captions 



Fig.l Renormalization group flows for the cubic model in three dimensions. 

Fig. 2 Log(Dn) as a function of Log(L). Triangles squares and circles denote the 
N = 4, N = 3 and N = 2 data respectively. Errors are not reported since they 
are smaller than the symbol sizes. To render easier the comparison among 
the three sets of data, all the values of Dr have been normalized to the best 
fit value of the constant C (see tab. 9 for the value of C). The three lines 
correspond to the best fits obtained neglecting the L = 6 derivative. 

Fig. 3 The N = 3 data only, plotted with a much higher resolution. The dotted line 
corresponds to the best fit including L = 6, while the dashed line corresponds 
to the L = [8 — 32] fit. All the points are normalized as in fig. 2. 
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